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A recently introduced phenomenological model to simulate flow control applications 
using plasma actuators has been further developed and improved in order to expand its use 
to complicated actuator geometries. The new modeling approach eliminates the requirement 
of an empirical charge density distribution shape by using the embedded electrode as a 
source for the charge density. The resulting model is validated against a flat plate 
experiment with quiescent environment. The modeling approach incorporates the effect of 
the plasma actuators on the external flow into Navier Stokes computations as a body force 
vector which is obtained as a product of the net charge density and the electric field. The 
model solves the Maxwell equation to obtain the electric field due to the applied AC voltage 
at the electrodes and an additional equation for the charge density distribution representing 
the plasma density. The new modeling approach solves the charge density equation in the 
computational domain assuming the embedded electrode as a source therefore automatically 
generating a charge density distribution on the surface exposed to the flow similar to that 
observed in the experiments without explicitly specifying an empirical distribution. The 
model is validated against a flat plate experiment with quiescent environment. 
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Nomenclature 

electric field, N/C 
permittivity, e — £ r £ a 
relative permittivity 

permittivity of free space, 8.854xl0" 12 C 2 /Nm 2 

body force vector, N/m 3 

total electric potential. Volt, <f> = <j> + cp 

electric potential due to external electric field. Volt 

electric potential due to net charge density. Volt 

length of the electrode 

Debye length, m 

location parameter for Gaussian distribution 
frequency, Hz 

Reynolds number based on inlet velocity and axial chord 
net charge density. C/m 3 
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t 

= time, s 

u in 

= inlet velocity 

u 

= streamwise velocity 

x, y 

= coordinates 

y n 

= wall normal distance 


I. Introduction 

P lasma actuators have been shown to be effective in several flow-control-related applications including flow 
separation and boundary layer control. The plasma actuators consist of two electrodes that are located on a 
surface separated by a dielectric material as shown in Fig. 1. A high-voltage AC supplied to the electrodes causes 
the air in their vicinity to weakly ionize. The ionized air (plasma), in the presence of the electric field gradient 
produced by the electrodes, result in a body force vector acting on the external flow that can induce steady or 
unsteady velocity components. The effectiveness of plasma actuators in controlling flow separation has been 
demonstrated by several researchers [1-9]. These experiments showed that a range of parameters have to be taken 
into consideration for effective flow control including the location of the actuators on the surface, orientation, size, 
and relative placement of the embedded and exposed electrodes, applied voltage, and frequency of the actuation. 
Due to a large number of parameters involved, optimizing the performance of actual applications can be a fairly 
complicated task. CFD simulations can provide a useful tool in design and optimization of such complex flow 
control systems. 

Computational studies of plasma flow control have been limited in comparison to the vast number of 
experimental studies reported [10-14]. Recently, the authors developed a simple phenomenological numerical 
simulation methodology for active flow separation control applications using plasma actuators [15]. In this new 
approach, the effect of the plasma actuators on the external flow is incorporated into Navier Stokes computations as 
a body force vector. The body force is obtained as a product of the net charge density and the electric field. The new 
model solves the Maxwell equation to obtain the electric field due to the applied AC voltage at the electrodes and an 
additional equation for the charge density representing the plasma density. However, there was an assumption 
made to the boundary condition of the charge density on the surface of the above the embedded electrode and this 
makes the extension of the model to complex and/or 3D flows impossible. This paper is an attempt to remedy this 
weakness. The model is summarized in the next section without the readers having to refer to our previous paper. 

II. Modeling Plasma Actuator Physics in CFD Computations 


Plasma 



Figure 1. Schematic of plasma actuator. 

The body force that the plasma actuator induces on the external flow can be expressed in terms of the applied 
voltage and incorporated into the Navier Stokes equations. By neglecting magnetic forces, the electrohydrodynamic 
(EHD) force can be expressed as 

fs=P c E ( 1 ) 

where, f B is the body force per unit volume, p c is net the charge density and E is the electric field. If the time 
variation of the magnetic field is negligible, as is often the case in plasma, the Maxwell’s equations give rise to 
Vx E ~ 0. This implies that the electric field can be derived from the gradient of a scalar potential [16]: 

E = -V<P. (2) 

Gauss’s law yields: 

V-(eE)=p c (3) 


2 

American Institute of Aeronautics and Astronautics 


or 


(4) 


V • (fV<5) = -p c 

where e is the permittivity. The permittivity can be expressed as: 

e = e r e 0 (5) 

where £ r i s the relative permittivity of the medium, and £ o is the permittivity of free space. 

If we write the net charge density within the plasma at any point in terms of the elementary charge e , the 
background plasma density n o , the temperature of the species T, and introduce the Debye length A d the charge 
density can be expressed in terms of the potential <f> and the Debye length A d [16], 

Pc /£ o =(- ( 6 ) 

The body force given by Eq. (1) can be calculated using Eqs. (2) and (6) by imposing the applied voltage at the 
electrodes as boundary conditions. 

In the current modeling approach the electric potential is decomposed into two parts by using the superposition 
technique [15]. Since the gas particles are weakly ionized, we can assume the potential <t> can be decoupled into two 
parts: one being a potential due to the external electric field, 0 , and the other being a potential due to the net charge 
density in the plasma, cp , 

<& = 0+(p. (7) 

This approach is similar to the one used in numerical simulation of electroosmotic flows in which case the 
external electric field generates a force on the charged particles creating flow[17]. 

If we assume that the Debye thickness is small and the charge on the wall is not large, the distribution of charged 
species in the domain is governed by the potential caused by the electric charge on the wall and is largely unaffected 
by the external electric field. Therefore, we can write two independent equations in terms of these two potentials, 
one for the external electric field due to the applied voltage at the electrodes: 


V ■ (eV 0>) = 0 

(8) 

and another one for the potential due to the charged particles: 


V(£ r V (p)=-(p c l £ o ). 

(9) 

Using pj£ 0 =(— 1/AjW Eq. (9) can be written as[16]. 


V {eyp c ) = Pj 

(10) 


Eq. (8) provides the solution for the electric potential, 0 , using the applied voltage on the electrodes as boundary 

conditions and p ( is obtained from Eq. (10) with a prescribed boundary condition on the surface over the 

embedded electrode. The boundary condition is synchronized with the applied voltage on the electrode and the 
resultant body force vector is computed by 

f B =p c E = p c (-V0). (11) 

The boundary conditions and the computational domains of Eqs. (8) and (10) for a single pair of electrodes are 
shown in Figures 2a and 2b, respectively. Equation (8) is solved for the electric potential, 0 , imposing the applied 

voltage on the electrodes as boundary conditions and using the appropriate £ r value on both the air side and the 
wall. For air, £ r =1.0 and the experiments we considered for test cases used Kapton as the dielectric material which 
has an £ r value of 2.7. On the wall-air interface harmonic mean of £ rl and £ rl must be used in order to conserve 
electric field[ 15]. 
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Air Side: 

On Outer Boundaries: 

Air Side: 

On Outer Boundaries: 

Vie rl vp) = 0 

^ = 0 

V(£yp c ) = p c iX- d 

o 

II 


dn 

£,i = 1.0 


ii 

£ r = mean(£ rl ,£ rl ) 

A d = 0.00017m 

£ r =mean(£ rl ,£ r2 ) 


Dielectric Material: 

' 

V(£ r2 vp) = 0 
£, 2 = 2.7 

p = 0 


V (£, V A) = /V^ C 


£ r2 = 2.1 


Pc=Pc(t) 


(a)Boundary conditions for Eq. (8) 


(b) Boundary conditions for Eq. (10) 


Figure 2. Boundary conditions and computational domain for Eqs. (8) and (10). 


Applied AC voltage imposed at the exposed (upper) electrode as boundary condition is: 

(i2) 

The wave form function, /(f), can be a sine wave given by: 

fit) = sin(2^<yf) (13) 

where co is the frequency and (f . ) max is the amplitude, both of which are known quantities from experiment. The 
embedded electrode is prescribed as ground by setting the electric potential to zero on that electrode. At the outer 
boundaries, dp/dn = 0 is assumed as shown in Figure 2a. 

In order to obtain a net charge density distribution over the embedded electrode resembling the plasma distribution 
observed in experiments [16,18,19], we introduce a new approach which eliminates the need for prescribing a 
charge density distribution on the surface over the embedded electrode. In this new method, Eq. (10) is solved 
assuming the lower electrode as a source for the net charge density, p c , which is synchronized with the time 
variation of the applied voltage pit) on the exposed electrode, /(f) in equation (12): 

p c it)=prm (i4) 

where pf m is the maximum value of the charge density allowed in the domain (in Coulomb/m 3 ) and is a parameter 
to be determined later. By assuming the embedded electrode as a source for the charge density, solution of Eq. (10) 
automatically results in a charge density distribution on the surface above the embedded electrode similar to the 
distribution observed in experiments and the model does not require prescription of a charge density distribution on 
the surface and thus makes the modeling of generic actuator geometries including generic shapes and orientations of 
electrodes straight forward. 

At this stage, it should be noted that in order to solve Eq. (10) we need to specify two parameters, namely, 
p'f" and A d . These parameters control the strength of the plasma actuator’s effects on the flowfield and extent of 
these effects into the flowfield. The Debye length, A d , is in the order of 0.00017m for the case considered here. The 
maximum charge density, p" r ‘" , is a parameter to be calibrated using available experimental data. 

It is important to note that, Eqs. (8) and (10) need to be solved only once at the beginning of Navier Stokes 
computations since these equations do not contain a time derivative term. This can be achieved by writing these 
equations in nondimensional forms as illustrated in Figs. 3a and 3b. The quantity p can be normalized by the value 

of AC voltage of the exposed electrode, pit), in equation (12). Equation (8) can then be solved by imposing a 
constant boundary condition equal to unity at the upper electrode. Once the dimensionless p distribution is 
determined, the dimensional p values at any given time can be obtained by multiplying this distribution with the 
corresponding value of pit ) given by equation (12). Similarly, Eq. (10) can be solved only once at the beginning of 
computations by using the dimensionless equation for the charge density distribution, which is normalized by 
Pf" fit) ■ This implies that the boundary condition for the dimensionless charge density on the embedded electrode 
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is equal to unity. Once the solution for the dimensionless charge density is established, the dimensional values at any 
time can be calculated from this distribution by multiplying it with the corresponding value of p™ x f(t) . 


This modeling approach for plasma actuators is implemented in the GHOST code developed at University of 
Kentucky. GHOST is a pressure-based code based on SIMPLE algorithm with second order accuracy in both time 
and space. This code is capable of handling complex geometries, moving and overset grids, and includes 
multiprocessor computation capability using MPI. The overset grid capability of the code enables incorporation of 




V.(fV f) = 0 V • (e r Vp* c ) = pj A 2 d 


A (0 = A max /(0 
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(a) Dimensional variables (b) Nondimensional forms of Eqs. (8) and (10) 

Figure 3. Dimensional and nondimensional variables and nondimensional form for Eqs. (8) and (10). 

plasma actuators into the computations with relative ease since electrodes can be defined as individual solid blocks. 
The domain can be divided into two separate computational domains: one for the air side and the other for the 
dielectric wall. The GHOST code has been previously validated against a wide range of test cases and flow 
conditions and has been used extensively in several low pressure turbine related publications[20,21]. 

III. Model Calibration and Validation Using Quiescent Flow Experiments 

In order to calibrate the maximum charge density on the wall, /T"” , appearing in the model we employed the 
quiescent flow experiments (Ref. 22) conducted using a single pair of electrodes to better isolate the effects of the 
actuator on the surrounding air. The details of the actuator geometry and experimental set up are given in Reference 
[15]. The actuator consists of two 10mm wide, 0.102mm thick conductive copper strips as electrodes which are 
separated by a 0.127mm thick Kapton dielectric with a£ r value of 2.7. Streamwise spacing of electrodes is 0.5mm. 

In the experiments the lower electrode was grounded and plasma region was generated using a square wave with 
frequency of, 0 ) = 4.5 kHz and amplitude of =5 kV. It should be noted that the experimental data is 
preliminary and it is used here only to demonstrate the proof of concept for the modeling approach. In the 
computations our aim was to match the maximum velocity observed in the experiments as well as the 
experimentally observed flow pattern shown in Figure 4a. From the experiment it was observed that the flow was 
drawn into the surface region above the embedded electrode by the plasma induced body force. This resulted in a jet 
issuing to the right of the actuator with a maximum velocity of approximately 1 m/s. 

Based on the flow pattern and maximum velocity criteria, p"'“ K =0.00750 C/m 3 seems to agree well with the 
experiment. The streamlines obtained from the computation using these values are shown in Figure (4b). The 
computed flow pattern compares favorably with the experimental flowfield shown in Figure (4a). The computed 
maximum velocity is lm/s also matching the experimental value. It should be noted that the boundary layer 
obtained by the computation appears to be thinner than the experimental data. 
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x(cm) 


(a) Experiment[22] (b) Computation 

Figure 4. Comparison of experimental[22] and computed streamlines for the plasma actuator in 
quiescent flow (The actuator interface is located at the 2cm tick.) 


The computed electric potential distribution in the vicinity of the electrodes obtained from equation (8) is shown 
in Figure 5 along with the streamlines of the actuator induced flow. The computed electric potentials show that the 
strongest electric potential variation, or the electric field, is in the region between the two electrodes. This is also the 
region where the strongest concentration of plasma is observed in the experiments. The streamlines indicate that the 
flow is pulled from above into this region and jetted to the right direction as observed experimentally. 



Figure 5. Computed electric potential contours and streamlines in the vicinity of the electrodes 

The computed charge density distribution obtained from solution of Eq. (10) is shown in Figure 7 on the air 
side. It must be noted that this charge density distribution over the surface is obtained without any explicit 
specification of a charge density variation. Using the embedded electrode as a source for Eq. (10) results 
automatically in the charge density distribution shown in Figure 7, eliminating any need for explicit specification of 
a charge density distribution. This aspect of the model is important since it makes the model applicable to any 
generic electrode shape or orientation without the need of explicit prescription of a charge density distribution. 

The magnitude of the resulting body force computed from Eq. (11) is shown in Figure 8. The maximum body 
force is concentrated in the region over the upper left corner of the embedded electrode where the electric field is the 
strongest. This is consistent with the experimental observations. 
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Figure 7. Computed charge density contours in the 
vicinity of the electrodes 


Figure 8. Computed body force magnitude contours 
in the vicinity of the electrodes 


These results show that the modeling technique developed can reproduce the effects of plasma actuators 
observed in experiments and the approach is promising in computation of plasma flow control applications 
involving generic actuator geometries and orientations. 

IV. Concluding Remarks 

A recently developed phenomenological model to simulate flow control applications using plasma actuators has 
been further developed and improved in order to expand its applicability to complicated actuator geometries. The 
new model eliminates the requirement of an empirical charge density distribution shape by using the embedded 
electrode as a source for the charge density. The resulting model is validated against a flat plate experiment with 
quiescent environment. The modeling approach incorporates the effect of the plasma actuators on the external flow 
into Navier Stokes computations as a body force vector which is obtained as a product of the net charge density and 
the electric field. The model solves the Maxwell equation to obtain the electric field due to the applied AC voltage at 
the electrodes and an additional equation for the charge density distribution representing the plasma density. These 
two equations need to be solved only once before the computation. The new modeling approach solves the charge 
density equation in the computational domain assuming the embedded electrode as a source therefore automatically 
generating a charge density distribution over the embedded electrode. This makes the model applicable to any 
generic electrode shape or orientation without the need of explicit prescription of a charge density distribution. 
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